Proyecto final, parte 1

Movimiento browniano geométrico

Valoración de opciones para el precio de Micron Techonology

Juliana Vallejo, Paola Fernández, María Camila Vásquez

In [1]:
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
from scipy import stats
from statsmodels.stats.diagnostic import lilliefors
import warnings
warnings.filterwarnings("ignore")
%matplotlib inline
import plotly.offline as py
import plotly.graph_objs as go
py.init_notebook_mode(connected=True)

Lectura de datos

In [2]:
X = pd.read_csv('datos.csv')['MU US Equity']
X = X.reindex(index=X.index[::-1])
X = np.array(X)
N = len(X)
R = np.zeros(N)
# Retornos instantáneos
for j in range(1,N):
        R[j-1] = (X[j] - X[j-1])/X[j-1]

Gráfica de la acción

In [3]:
trace1 = go.Scatter(
    y = X,
    x = np.arange(N),
    name = '$X_1$'
)

data = [trace1]
layout = go.Layout(
    title=go.layout.Title(
        text='Precio de Micron Techonology',
        xref='paper',
        x=0
    ),
    xaxis=go.layout.XAxis(
        title=go.layout.xaxis.Title(
            text='Tiempo',
            font=dict(
                family='Arial',
                size=18,
                color='#7f7f7f'
            )
        )
    ),
    yaxis=go.layout.YAxis(
        title=go.layout.yaxis.Title(
            text='$X_t$',
            font=dict(
                family='Arial',
                size=18,
                color='#7f7f7f'
            )
        )
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='Accion')

Gráfica de los retornos instantáneos

In [4]:
trace1 = go.Scatter(
    y = R,
    x = np.arange(N),
    name = '$R_t$'
)

data = [trace1]
layout = go.Layout(
    title=go.layout.Title(
        text='Retornos instantáneos',
        xref='paper',
        x=0
    ),
    xaxis=go.layout.XAxis(
        title=go.layout.xaxis.Title(
            text='Tiempo',
            font=dict(
                family='Arial',
                size=18,
                color='#7f7f7f'
            )
        )
    ),
    yaxis=go.layout.YAxis(
        title=go.layout.yaxis.Title(
            text='$R_t$',
            font=dict(
                family='Arial',
                size=18,
                color='#7f7f7f'
            )
        )
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='Retornos')
In [5]:
trace1 = go.Histogram(
    x = R
)

data = [trace1]
layout = go.Layout(
    title=go.layout.Title(
        text='Histograma Retornos instantáneos',
        xref='paper',
        x=0
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='Retornos')

Test de normalidad

In [6]:
jb = stats.jarque_bera(R)
shapiro = stats.shapiro(R)
In [7]:
df = pd.DataFrame()
df['Jarque Bera'] = [jb[1]]
df['Shapiro Wilks'] = [shapiro[1]]
In [8]:
def hFD(a, k_max): 
# Higuchi FD
    L = []
    x = []
    N = len(a)

    for k in range(1,k_max):
        Lk = 0
        for m in range(0,k):
            #we pregenerate all idxs
            idxs = np.arange(1,int(np.floor((N-m)/k)),dtype=np.int32)
            Lmk = np.sum(np.abs(a[m+idxs*k] - a[m+k*(idxs-1)]))
            Lmk = (Lmk*(N - 1)/(((N - m)/ k)* k)) / k
            Lk += Lmk

        L.append(np.log(Lk/(m+1)))
        x.append([np.log(1.0/ k), 1])

    (p, r1, r2, s)=np.linalg.lstsq(x, L)
    return p[0]
In [9]:
def hurst(signal):
    tau = []; lagvec = []

    #  Step through the different lags
    for lag in range(2,20):

    #  produce price difference with lag
        pp = np.subtract(signal[lag:],signal[:-lag])

    #  Write the different lags into a vector
        lagvec.append(lag)

    #  Calculate the variance of the difference vector
        tau.append(np.std(pp))

    #  linear fit to double-log graph (gives power)
    m = np.polyfit(np.log10(lagvec),np.log10(tau),1)

    # calculate hurst
    hurst = m[0]

    return hurst
In [10]:
df['Coeficiente de Hurst'] = hurst(X)
df['Dimensión Fractal'] = hFD(X, 8)
In [11]:
print(df.to_latex(index=False))
\begin{tabular}{rrrr}
\toprule
 Jarque Bera &  Shapiro Wilks &  Coeficiente de Hurst &  Dimensión Fractal \\
\midrule
    0.066684 &       0.085535 &              0.429953 &           1.494463 \\
\bottomrule
\end{tabular}

In [12]:
from statsmodels.graphics.tsaplots import plot_pacf
pacf = plot_pacf(X,lags=10)
plt.title("PACF para Micron Techonology")
pacf.show()

Estimación de parámetros

In [13]:
dt = 1/len(R)
media = np.mean(R)
varianza = np.var(R)
mu = media/dt
sigma = np.sqrt(varianza/dt)
print('$\mu: '+str(mu)+'$\n'+'$\sigma: '+str(sigma)+'$')
$\mu: -0.048218673042336915$
$\sigma: 0.4831365452648228$
In [14]:
ms = []
In [15]:
# Parametros
dt = 1/252
#mu = 0.5
sg = sigma
# Condicion inicial
X0 = X[0]
k = 1000
S = np.zeros((k,N))
for i in range(k):
    S[i,0] = X0
for i in range(k):
    errr = []
    for j in range(1,N):
        S[i,j] = X[j-1] + mu*X[j-1]*dt + sg*X[j-1]*np.sqrt(dt)*np.random.normal()
        errr.append(np.abs(S[i,j]-X[j])/X[j])
    ms.append(np.mean(errr))
In [16]:
print(np.mean(ms))
0.034146221938439004
In [17]:
data=[]
for i in range(k):
    trace1 = go.Scatter(
        y = S[i],
        x = np.arange(N),
        showlegend=False
    )
    data.append(trace1)

trace2 = go.Scatter(
    y = X,
    x = np.arange(N),
    mode = 'lines',
    line = dict(color = 'rgb(0,0,0)'),
    name = 'Real'
)
data.append(trace2)
layout = go.Layout(
    title=go.layout.Title(
        text='Simulación con parámetros estimados',
        xref='paper',
        x=0
    ),
    xaxis=go.layout.XAxis(
        title=go.layout.xaxis.Title(
            text='Tiempo',
            font=dict(
                family='Arial',
                size=18,
                color='#7f7f7f'
            )
        )
    ),
    yaxis=go.layout.YAxis(
        title=go.layout.yaxis.Title(
            text='$X_t$',
            font=dict(
                family='Arial',
                size=18,
                color='#7f7f7f'
            )
        )
    )
)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, filename='Estimacion')

Opciones

In [18]:
# Parametros
# condicion inicial
S0 = X[-1]
# Tasa de interés libre de riesgo
r = 0.05
# Precio de ejercicio
K = 1.1*S0
In [19]:
op = pd.DataFrame()
In [20]:
def black(T):
    d1 = (np.log(S0/K) + (r+0.5*sigma**2)*T)/(sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    import scipy.stats as ss
    N = ss.norm(0,1)
    fcall = S0*N.cdf(d1) - K*np.exp(-r*T)*N.cdf(d2)
    return fcall    
In [21]:
# Simulación montecarlo
def montecarlo(T):
    k = 100
    M = 200
    N = 250
    dt = T/N
    S = np.zeros([k,N]) # Filas cada tracyectoria, columnas el tiempo
    payoff=np.zeros(k) #en cada simulacion saco el payoff de cada trayectoria
    opcion=np.zeros(M) #Tantas opciones como simulaciones
    for h in range(M):
        for i in range(k):
            S[i,0] = S0
        for i in range(k):
            for t in range(N-1):
                S[i,t+1]=S[i,t]+r*S[i,t]*dt+sigma*S[i,t]*np.sqrt(dt)*np.random.normal()
            payoff[i]=np.max([S[i,N-1]-K,0]) 
        opcion[h]=np.exp(-r*T)*np.mean(payoff); 
    opcionmedio = np.mean(opcion)
    opcionmax = np.max(opcion)
    opcionmin = np.min(opcion)
    return opcionmedio
In [22]:
# Árboles binomiales
def binomial(T):
    n = 50
    dt = T/n
    u = np.exp(sigma*np.sqrt(dt)) # Proporción de subida
    d = 1/u # Proporción de bajada
    p = (np.exp(r*dt)-d)/(u-d) # Probabilidad de salto
    f = np.zeros((n+1,n+1))
    for j in range(n+1):
        a = np.maximum(S0*(u**(n-j))*(d**(j)) - K,0)
        f[n,j] = a
    for i in range(n-1,-1,-1):
        for j in range(i+1):
            f[i,j] = np.exp(-r*dt)*((1-p)*f[i+1,j+1] + (p)*f[i+1,j])
    return f[0,0]
In [23]:
# Diferen
def diferencias(T):
    S_max = 2*S0 
    dS = 0.2   
    M = int(np.ceil(S_max//dS))
    dT_tmp = (dS/(sigma*S_max))**2
    N_tmp = T/dT_tmp
    N = int(np.ceil(N_tmp))
    dT = T/N
    a = np.zeros(M+1)
    b = np.zeros(M+1)
    c = np.zeros(M+1)
    for j in range(M+1):
        a[j] = (dT/(1+r*dT)) * ((sigma**2 * j**2)/2 - (r * j/2))
        b[j] = (dT/(1+r*dT)) * ((1/dT) - sigma**2 * j**2)
        c[j] = (dT/(1+r*dT)) * ((sigma**2 * j**2)/2 + (r * j/2))
    # Crear Malla F
    F = np.zeros((N+1, M+1))

    # Condiciones de frontera
    F[:,0] = 0

    for j in range(M):
        F[N, j] = max(j*dS - K, 0)

    F[:, M] = S_max - K

    for i in reversed(range(N)):
        for j in reversed(range(1, M)):
            F[i,j] = a[j] * F[i+1, j-1] + b[j] * F[i+1, j] + c[j] * F[i+1, j+1]
    P = S0//dS;
    FF = F[0,int(P)]
    return FF
In [24]:
aa = [1/4, 1/2, 1, 1/6]
for t in aa:
    rrr = []
    rrr.append(black(t))
    rrr.append(montecarlo(t))
    rrr.append(binomial(t))
    rrr.append(diferencias(t))
    rrr = pd.Series(rrr)
    op[t] = rrr
In [25]:
op.index = ['Black Scholes', 'Simulación Montecarlo', 'Árboles binomiales', 'Diferencias Finitas']
In [26]:
print(op.to_latex())
\begin{tabular}{lrrrr}
\toprule
{} &  0.250000 &  0.500000 &  1.000000 &  0.166667 \\
\midrule
Black Scholes         &  2.668580 &  4.541519 &  7.335854 &  1.888105 \\
Simulación Montecarlo &  2.674308 &  4.458363 &  7.172518 &  1.857926 \\
Árboles binomiales    &  2.685097 &  4.514264 &  7.361655 &  1.897326 \\
Diferencias Finitas   &  2.664212 &  4.527170 &  7.211015 &  1.884434 \\
\bottomrule
\end{tabular}

In [ ]: